#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

          //////////////////
          // Input   ///////
          /////////////////

	  int NOa2;
	  int Nf;
          double beta, Uplaq;
          char dir[1024];
          char gauge_fixing[1024];
          char gauge_action[1024];



          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-Nf",3)==0)  sscanf(argv[i],"-Nf=%d",&Nf);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-gauge_fixing",13)==0)  sscanf(argv[i],"-gauge_fixing=%s",gauge_fixing);
            if(strncmp(argv[i],"-gauge_action",13)==0)  sscanf(argv[i],"-gauge_action=%s",gauge_action);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta",5)==0)  sscanf(argv[i],"-beta=%lf",&beta);
            if(strncmp(argv[i],"-Uplaq",6)==0)  sscanf(argv[i],"-Uplaq=%lf",&Uplaq);
            if(strncmp(argv[i],"-NOa2",5)==0)  sscanf(argv[i],"-NOa2=%d",&NOa2);

          }

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  Equivalent_Momenta(SinP2);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(nn);

  // IO

   char file[1024];
   char Oa2print[1024];
   char dir_qav[1024];
   char output[1024];
   sprintf(dir_qav,"%s/QAV/",dir);

    //////////////////////////////////
    // Read Zq 
     /////////////////////////////

    double ***Zq=malloc(sizeof(double**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) Zq[imom]=malloc(sizeof(double*)*(NOa2));
    for (int imom=0; imom<Nmom_eq; imom++) for(int iOa2=0; iOa2<NOa2; iOa2++) Zq[imom][iOa2]=malloc(sizeof(double)*(Nclust+1));

    O55 ***DpPC=malloc(sizeof(O55**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) DpPC[imom]=malloc(sizeof(O55*)*(NOa2));
    for (int imom=0; imom<Nmom_eq; imom++) for(int iOa2=0; iOa2<NOa2; iOa2++) DpPC[imom][iOa2]=malloc(sizeof(O55)*(Nclust+1));

    double *aux_clust=malloc(sizeof(double)*(Nclust+1));

    sprintf(file,"Zq");
    for (int imom=0; imom<Nmom_eq; imom++) {
	    read_input_onemomenta(dir_qav,file,aux_clust,imom);
	    for (int iclust=0; iclust<Nclust+1; iclust++) Zq[imom][0][iclust]=aux_clust[iclust];
     }

    for (int imom=0; imom<Nmom_eq; imom++) {
	for (int ioperator=0; ioperator<5; ioperator++) for (int iprojector=0; iprojector<5; iprojector++)
          {
		sprintf(file,"DpPC%d%d",ioperator+1,iprojector+1);
		read_input_onemomenta(dir_qav,file,aux_clust,imom);
		for (int iclust=0; iclust<Nclust+1; iclust++) DpPC[imom][0][iclust][ioperator][iprojector]=aux_clust[iclust];
	  }
    }    



     ///////////////////////////////
     // Compute O(a2g2) corrections
     ///////////////////////////////


     //////////////// substracting terms     
     double **CorOa2_Zq=malloc(sizeof(double*)*(Nmom_eq));
     for (int imom=0; imom<Nmom_eq; imom++) CorOa2_Zq[imom]=malloc(sizeof(double)*(NOa2));

     O55 **CorOa2_D=malloc(sizeof(O55*)*(Nmom_eq));
     for (int imom=0; imom<Nmom_eq; imom++) CorOa2_D[imom]=malloc(sizeof(O55)*(NOa2));


	double  b0;    
    for (int imom=0; imom<Nmom_eq; imom++) 
    {
     //iOa2=0 -> uncorrected
     CorOa2_Zq[imom][0]=0.;
     for (int ioperator=0; ioperator<5; ioperator++) for (int iprojector=0; iprojector<5; iprojector++)  CorOa2_D[imom][0][ioperator][iprojector]=0.;
     //iOa2=1 corrected
     //iOa2=2 corrected boosted
     double g2;
     for (int iOa2=1; iOa2<NOa2; iOa2++)
     {
	  // Physical Review D Volume 48, Number 5, Viability of lattice perturbation theory, Lepage and Mackenzie
          if(iOa2==1) g2=6.0/beta;
          if(iOa2==2) g2=(6.0/beta)/(Uplaq);
	  // if(iOa2==2) printf("%8.5f\n ",g2);
	  if(iOa2==3) { g2=(6.0/beta)/(Uplaq) ; b0=(11*Nc-2*Nf)/3;  g2=pow(1/g2 -b0/(16*PI*PI)*pow(log(PI/sqrt(SinP2_eq[imom])),2)-4.702/(4*PI),-1);}
	   //if(iOa2==3) printf("%8.5f\n",g2);
          CorOa2_Zq[imom][iOa2]=correction_Zq_Oa2g2(g2,SinP2_eq[imom],SinP4_eq[imom],nn[imom],gauge_action,gauge_fixing);
	  correction_D_Oa2g2(CorOa2_D[imom][iOa2],g2,SinP2_eq[imom],SinP4_eq[imom],nn[imom],gauge_action,gauge_fixing);

     }
    }

    ////////////////// Zq corrected
    for (int imom=0; imom<Nmom_eq; imom++) 
       for (int iclust=0; iclust<Nclust+1; iclust++) 
         for (int iOa2=1; iOa2<NOa2; iOa2++) 
		Zq[imom][iOa2][iclust]=Zq[imom][0][iclust]-CorOa2_Zq[imom][iOa2];

    ////////////////// D corrected
    // Note that following correction_4f.nb the first index  in the corrector parameter corresponds to the projector and the second to the operator
    // that is: [Projector, LGreen]
    for (int imom=0; imom<Nmom_eq; imom++)
       for (int iclust=0; iclust<Nclust+1; iclust++)
         for (int iOa2=1; iOa2<NOa2; iOa2++)
 	   for (int ioperator=0; ioperator<5; ioperator++) for (int iprojector=0; iprojector<5; iprojector++) 
                DpPC[imom][iOa2][iclust][ioperator][iprojector]=DpPC[imom][0][iclust][ioperator][iprojector]-CorOa2_D[imom][iOa2][iprojector][ioperator];


     ///////////////////////
     /// WRITING OUTPUT
     ///////////////////////

      //printf(" WRITING OUTPUT IN %s\n",dir);
     double **aux_allmomenta_real=malloc(sizeof(double*)*(Nmom_eq));
     for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta_real[imom]=malloc(sizeof(double)*(Nclust+1));

   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
	 whichOa2(iOa2,Oa2print);
         sprintf(output,"%s/Z_ALLP/Zq_%s.out",dir_qav,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=Zq[imom][iOa2][iclust];
         write_output_allmomenta_eq_real(output,aux_allmomenta_real,SinP2_eq);

         sprintf(file,"Zq_%s",Oa2print);
         write_output_momentabymomenta_eq_real(dir_qav,file,aux_allmomenta_real,SinP2_eq);
   }

   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
        whichOa2(iOa2,Oa2print);
	for (int ioperator=0; ioperator<5; ioperator++) for (int iprojector=0; iprojector<5; iprojector++)
	{
         sprintf(output,"%s/D_ALLP/DpPC%d%d_%s.out",dir_qav,ioperator+1,iprojector+1,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=DpPC[imom][iOa2][iclust][ioperator][iprojector];
         write_output_allmomenta_eq_real(output,aux_allmomenta_real,SinP2_eq);

         sprintf(file,"DpPC%d%d_%s",ioperator+1,iprojector+1,Oa2print);
         write_output_momentabymomenta_eq_real(dir_qav,file,aux_allmomenta_real,SinP2_eq);
	}
   }


    free(Zq); free(DpPC);
    return 0;

}

